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Abstract. A third order shock-capturing numerical scheme for three-dimensional special relativistic magnetohy- 
drodynamics (3-D RMHD) is presented and validated against several numerical tests. The simple and efficient 
central scheme described in Paper I (Del Zanna and Bucciantini, Astron. Astrophys., 390, 1177-1186, 2002) for 
relativistic hydrodynamics is here extended to the magnetic case by following the strategies prescribed for classical 
MHD by Londrillo and Del Zanna (Astrophys. J., 530, 508-524, 2000). The scheme avoids completely spectral 
decomposition into characteristic waves, computationally expensive and subject to many degenerate cases in the 
magnetic case, while it makes use of a two-speed Riemann solver that just require the knowledge of the two local 
fast magnetosonic velocities. Moreover, the onset of spurious magnetic monopoles, which is a typical problem 
for multi-dimensional MHD upwind codes, is prevented by properly taking into account the solenoidal constraint 
and the specific antisymmetric nature of the induction equation. Finally, the extension to generalized orthogonal 
curvilinear coordinate systems is included, thus the scheme is ready to incorporate general relativistic (GRMHD) 
effects. 

Key words, magnetohydrodynamics (MHD) - relativity - shock waves - methods: numerical 



1. Introduction 

Most of the astrophysical sources of high-energy radiation 
and particles are believed to involve the presence of rela- 
tivistic motions in magnetized plasmas. For example, the 
radio emission from extra galactic jets (especially from ter- 
minal radio lobes) or from plerion-like supernova remnants 
is due to synchrotron radiation produced by relativistic 
electrons spiraling around magnetic field lines, thus indi- 
cating the presence of significant magnetic fields. Strong 
magnetic fields are supposed to play an essential role in 
converting the energy of accreting material around super- 
massive black holes at the center of Active Galactic Nuclei 
(AGNs), into powerful relativistic jets escaping along open 
field lines (Begelman et al. 1984). Similar phenomena may 
be at work in the galactic compact X-ray sources known 
as microquasars (Mirabel & Rodriguez 1994). These pro- 
cesses involve the interaction of relativistic gasdynamic 
flows and shocks with strong magnetic fields, which have 
now started to be studied via computer simulations (see 
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Meier et al. 2001 for a review). Powerful relativistic blast 
shocks should also be at the origin of the still mysteri- 
ous gamma-ray bursts (GRBs; Meszaros & Rees 1992). 
Moreover, the presence of magnetic field has been invoked 
in various astrophysical objects to explain both their mor- 
phology and evolution by applying simplified analytical 
models to basic plasma physics effects (e.g. the magnetic 
pinch and kink instabilities may affect the structure of 
both AGN and microquasar jets, and possibly also the 
overall shape of pulsar wind nebulae) , although a detailed 
study of the nonlinear and turbulent regimes is still lack- 
ing. 

Due to the extreme complexity and richness of the pos- 
sible effects arising in relativistic plasma physics, there is 
a very strong interest among the astrophysical commu- 
nity in the development of computer codes for both rela- 
tivistic hydrodynamics (RHD) and magnetohydrodynam- 
ics (RMHD), since in most cases only numerical simula- 
tions are able to cope with the evolution of such phenom- 
ena. After some early attempts based on non-conservative 
schemes that handled shocks with the aid of large arti- 
ficial viscosity and resistivity, it is only during the last 
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decade that conservative shock-capturing Godunov-type 
numerical codes, already successfully applied to gasdy- 
namic problems, have started to be applied to RHD too, 
achieving both high accuracy in smooth regions of the sim- 
ulated flow and sharp discontinuous profiles at shocks (e.g. 
Marquina et al. 1992; Schneider et al. 1993; Balsara 1994; 
Duncan & Hughes 1994; Eulderink & Mellema 1994; Font 
et al. 1994; Dolezal & Wong 1995; Falle & Komissarov 
1996; Donat et al. 1998; Aloy et al. 1999; Del Zanna & 
Bucciantini 2002). 

However, in spite of the success of Godunov-type RHD 
codes and, at the same time, of the presence of various 
extensions of gasdynamic schemes to classical MHD (see 
the recent review by Toth, 2000), to date only a cou- 
ple of RMHD schemes have been described in the litera- 
ture. Both codes are second order accurate and are based 
on linearized Rlemann solvers (Roe matrix) in the def- 
inition of fluxes at cell interfaces. This process involves 
the decomposition of primitive variables in a set of char- 
acteristic waves, each of them propagating a single dis- 
continuity, and a further composition to obtain the nu- 
merical upwind fluxes. Moreover, a certain amount of 
extra artificial viscosity is often needed to stabilize the 
schemes in particular degenerate cases. The two codes 
are described in Komissarov (1999a; KO from now on), 
which is a truly multidimensional scheme, and in Balsara 
(2001; BA from now on), the latter tested just against 
one-dimensional (1-D) shock-tube problems. There is actu- 
ally another RMHD code, which has been extensively used 
in relativistic 2-D and 3-D jet simulations (e.g. Koide et 
al. 1996; Nishikawa et al. 1998), later extended to general 
relativistic (GRMHD) effects (with given Schwarzschild or 
Kerr metrics) and applied to the jet formation mechanism 
(e.g. Koide et al. 1999; 2000). However, this code cannot 
be regarded as belonging to the Godunov family, since it is 
based on a second order Lax-Wcndroff scheme, thus with 
a very high level of implicit numerical viscosity. Moreover, 
a complete set of the standard numerical tests, needed to 
check the properties of any shock-capturing scheme, has 
never been published for such code, so it is difficult to com- 
ment on their results and to compare the respective code 
performances (especially on contact discontinuities, where 
shock-capturing codes not based on linearized Ricmann 
solvers are usually less accurate). 

The reasons behind the difficulty of extending shock- 
capturing relativistic gasdynamic codes to the magnetic 
case are essentially the same encountered in building clas- 
sical MHD schemes, but amplified, so to say, by the special 
relativistic effects. These difficulties may be summarized 
basically in two classes of problems. The first is concerned 
with the eigenstructure of the 1-D MHD system, which 
is much more complex than in the fluid case since now 
seven characteristic waves are involved and many different 
degeneracies may occur (depending on the relative orien- 
tation of the magnetic and velocity vectors). These prob- 
lems of non-strict hyperbolicity can be cured by accurate 
re- normalizations of the variables (Brio & Wu 1988; Roe 
& Balsara 1996), to assure their linear independence, and 



by introducing additional numerical viscosity in the de- 
generate cases. The second aspect is more crucial. The 
multidimensional MHD system, in conservative form, has 
a specific irreducible structure: the magnetic field (which 
is a pseudo-vector) is advanced in time by an antisymmet- 
ric differential operator, a curl, while all other variables are 
scalars or vectors advanced in time by a differential opera- 
tor of divergence form. We notice that this basic duality in 
the conservation laws of the MHD system is also fully in- 
variant under relativistic coordinate transformations: the 
covariant evolution equation for B splits into the classical 
induction equation and in the non-evolutionary solcnoidal 
V • B = condition. 

In numerical schemes where Godunov-type procedures, 
based on the divergence conservation form and on cell- 
centered variables, are also applied to the induction equa- 
tion, it comes out that magnetic field components develop 
unphysical discontinuities and numerical monopolcs which 
may grow in time. In RMHD this problem can be relevant: 
when the magnetic field is very strong, that is when the 
Alfven velocity approaches the speed of light, the various 
eigenvalues collapse one onto the other and it becomes 
very hard, from a numerical point of view, to distinguish 
among different physical states: thus, even very small er- 
rors in the definition of the magnetic field components, or 
in the flux derivatives (where the solcnoidal condition is 
implicitly assumed), often lead to unphysical states and 
to code crashing. Therefore, the proper character of the 
induction equation and the related preservation of the 
solenoidal constraint are fundamental issues in numerical 
RMHD. 

To cope with this class of numerical problems, two 
families of empirical solutions have been proposed: the 
cleaning methods, where the magnetic field components 
are re-defined at every time step (originally proposed by 
Brackbill & Barnes, 1980, it requires the solution of an 
additional Poisson equation), and the eight-wave method 
(Powell 1994), which modifies the MHD system by adding 
a new V • B variable, to be advected by the flow like 
the other quantities. These methods may in some cases 
alleviate (but not solve) the main difficulties. A more con- 
sistent way to handle this problem is given by the family 
of constrained transport (CT) methods (first introduced 
by Evans & Hawley, 1988), where the induction equation 
is correctly discretized to incorporate the solcnoidal con- 
straint as a main built-in property. Many schemes (Dai 
& Woodward 1998; Ryu et al. 1998; Balsara & Spicer 
1999b, to cite a few) take advantage of this method, but 
only in a restricted way, since all basic upwind procedures 
are still based on the standard (cell-centered) Godunov- 
type formalism and therefore the production of numerical 
monopoles is not avoided. A significant further advance 
has been proposed by KO: after an initial attempt to ex- 
tend the eight-wave method to RMHD (failed basically 
because of the above mentioned numerical problems), 
he finally turns to a CT scheme where the discretized 
divergence-free magnetic field components are correctly 
incorporated in the momentum-energy flux functions, al- 



L. Del Zanna et al.: A shock-capturing scheme for relativistic flows - II. 



though the upwind fluxes for the induction equations are 
still not properly defined, in our opinion. However, the 
only astrophysical application of such RMHD code pub- 
lished so far is the propagation of light relativistic jets 
embedded in a purely toroidal magnetic field (Komissarov 
1999b), where the solenoidal condition is actually auto- 
matically satisfied for simple geometrical reasons (the field 
is bound to remain always toroidal), thus this test is not 
stringent at all for the solenoidal condition preservation 
problem. 

To date, a fully consistent CT-based upwind scheme 
assuring exact V • B = condition has been proposed 
by Londrillo & Del Zanna (2000), LD for brevity. In fact, 
starting from a finite volume formulation of the solenoidal 
condition and of the induction equation, as in the original 
CT method, general recipes are given to reconstruct (to 
any order of accuracy) magnetic field variables at points 
needed for flux computations and to formulate approxi- 
mate Riemann solvers also for the induction equations in 
the CT form. Moreover, as an application, a third order 
Essentially Non-Oscillatory (ENO) central-type scheme 
was proposed and numerically validated against various 
tests. The so-called central schemes do not make use of 
time-consuming and system-dependent spectral decompo- 
sitions, and linearized Riemann solvers are replaced by 
Lax-Friedrichs type averages over the local Riemann fan. 
In this way, the only characteristic quantities entering the 
scheme are the local fastest velocities, and also the prob- 
lems related to the various degeneracies are thus avoided 
completely. The price to pay is just some additional nu- 
merical dissipation at contact and Alfvenic discontinuities, 
but the high order reconstruction is often able to compen- 
sate for these drawbacks. 

The central scheme described in LD was applied to the 
RHD system in Del Zanna & Bucciantini (2002), from now 
on simply Paper I, where the test simulations presented 
demonstrated the accuracy and stability of such scheme, 
even in highly relativistic situations, giving equivalent or 
better (thanks to its higher order) results than those pro- 
duced by much more elaborate Godunov-type algorithms. 
Here the same third order ENO-CT central scheme of LD 
is extended to the RMHD system, thus this paper may be 
considered as the generalization of Paper I to the magnetic 
case. Therefore, both the structure of the paper and the 
formalism used will be the same as in Paper I, to which 
the reader will be often referred, especially for some nu- 
merical scheme details or for test simulations comparisons. 
Finally, the CT scheme is extended to generalized orthog- 
onal curvilinear coordinates in the appendix, thus, the in- 
clusion of General Relativity effects with a given metric, 

1. e. the extension to GRMHD, may be easily achieved (see 
the appendix of Koide et al. 1999). 

2. Ideal RMHD equations 

The covariant fluid equations for special relativistic hy- 
drodynamics (RHD) were given in Paper I and here the 
same notation will be assumed throughout, that is all ve- 
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locities are normalized against the speed of light (c = 
1), Greek (Latin) indexes indicate four (three) vectors, 
g a P = diag{— 1, 1,1,1} is the Minkowski metric tensor (a 
flat space is assumed here for ease of presentation, for 
the extension to any set of orthogonal curvilinear coordi- 
nates see the appendix), and x a = (t,x J ) is the four vec- 
tor of space-time coordinates. The modifications needed 
to take electromagnetic forces into account are, like in 
classical MHD, the inclusion of extra terms in the energy- 
momentum conservation law and a new equation for the 
magnetic field, to be derived from Maxwell equations. Our 
derivation follows that of Anile (1989), also described in 
KO and BA. 

Written in terms of the (antisymmetric) electromag- 
netic tensor F af3 (F 0i = E i} F ij = B k with {i,j,k} = 
{1,2,3} and cyclic permutations) and of its dual F* a/3 — 
i e aff 7 5 F ^ s ( F *Qi = B h F* ij = -E k ), where e a ^ s is 
the Levi-Civita alternating pseudo-tensor, the covariant 
Maxwell equations are: 

3 a F a!i = -J f3 , 8 a F* afi = 0, (1) 

where J a is the four-current containing the source terms, 
constrained by the condition d a J a = 0, and we have as- 
sumed 47T — > 1. On the other hand, the electromagnetic 
contribution to the energy-momentum tensor is 

T»i = F a ~ f F l3 ~ l ~ ^F lS F^ s , (2) 

to be added to the fluid part in the conservation law 
d a T a P = 0. Finally, we must introduce the covariant rel- 
ativistic form of Ohm's law in the infinite conductivity 
approximation, E + v x B = 0, that translates into a 
condition of vanishing covariant electric field 

F af3 u = 0, (3) 

where u a = (7, 71/') is the four- velocity and 7 = u° = 
(1 — v 2 ) -1 / 2 is the Lorentz factor. Note that the other 
approximation needed to derive the classical MHD equa- 
tions, namely to neglect the displacement current, is not 
imposed in RMHD, of course, and the result is that the 
current, to be derived from the first of Eqs. ([!]), now 
depends on the time derivative of the electric field too: 
J = V x B d t E. 

The equations written so far are not easily compared 
with their MHD equivalent, due to the presence of the 
electromagnetic tensor and of its dual, both containing 
the electric field. However, thanks to Eq. (||), E may 
be substituted everywhere by defining a magnetic induc- 
tion four-vector as b a — F* af3 U/3, that allows to write 
the electromagnetic tensor in terms of u a and b a alone: 
F ~f5 = e <*p-y$ baU0 , The components of this new four- vector 
are 

b a = [ 7 (w-B),B/7 + 7(« -B)v], (4) 

and in the fluid comoving local rest frame we simply have 
b a = (0,B). Note the constraints u a b a — and \u\ 2 = 
u a u a = —1, so that \b\ 2 = b a b a > and b a is a space-like 
vector, with |6| 2 = B 2 /-/ 2 + (v ■ B) 2 . 
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Thanks to these definitions, the complete set of RMHD 
equations becomes: 

d a (pu a )=Q, (5) 

d a [{w + \b\ 2 )u a U P -b%P + {p+\b\ 2 /2)g al} ] = 0, (6) 

d a (u a b^-^b a ) = 0, (7) 

namely the equations of mass conservation, of total 
energy-momentum conservation, and of magnetic induc- 
tion. Here w = e + p is the relativistic enthalpy and 
e = p + p/(r — 1) is the relativistic energy per unit vol- 
ume for a T-law equation of state. Notice the analogy 
with classical MHD equations, easily obtained by letting 
v 2 <C 1, while RHD equations are recovered simply by 
letting b a = 0. 

2.1. Evolution equations and the V • B = constraint 

Godunov-type shock-capturing numerical methods devel- 
oped for classical Euler equations apply for any set of hy- 
perbolic conservation laws, and Paper I has shown the 
application to the RHD case. It is easy to verify that the 
equations for the fluid variables, Eqs. (||) to (||), retain the 
usual conservative form: 



du 
~dt 



2^, Q x i 

i=l 



o. 



(8) 



Here u is the vector of conserved variables and j l are 
their corresponding fluxes, along each direction, respec- 
tively given by 

u = [pu°, w tot u°u j - b°V,w tot u°u - b°b° - p tot ] T , (9) 

f = [pu\ wtotwV - VV + p tot 5V, w tot u°u l - b°b l f, (10) 

where we have defined Wtot = w+|b| 2 andptot — P+\b\ 2 

On the other hand, the equation for b a , Eq. ([?]), splits 
into two parts, which happen to be exactly the same as in 
classical MHD (this is not surprising since Maxwell equa- 
tions are Lorentz invariant). The spatial component gives 
the classical induction equation 

dB 



dt 



Vx£ = 0: 



E 



-v x B, 



(11) 



which is properly the time evolution equation for B. Note 
that the spatial differential operator is in a curl form, 
rather than in a divergence form as Eq. (JsJ) . This means 
that the evolution equation of each spatial component of 
B has a missing eigen-space, basically due to the anti- 
symmetry of the electromagnetic tensor in Eq. (Q) , as an- 
ticipated in the introduction. Thus, a total of just three 
independent magnetic fluxes (the electric field vector com- 
ponents, just one in 2-D) are needed for the evolution 
of B, while six independent fluxes were required for the 
momentum evolution. The other consequence of the ten- 
sor antisymmetric nature is that the time component of 
Eq. (@) becomes the usual MHD solenoidal constraint 



which is not an evolutionary equation but a differential 
constraint on the spatial derivatives of B. This constraint 
is usually regarded as just an initial condition, since the 
form of the induction equation assures its preservation in 
time. Therefore, also numerical schemes must be designed 
in a way that the specific divergence-free nature of the 
magnetic field is taken into account as a fundamental con- 
stitutive property, otherwise spurious magnetic monopoles 
will affect the overall solution and often the code stability 
itself. The CT schemes, and our specific implementation 
described in Sect. 3, are the class of numerical schemes 
based on this property. 

It is now apparent that Eqs. (|ll]) and (|l^) are substan- 
tially different from the evolutionary conservation laws 
in Eq. (|§|). This fundamental constitutional difference, in 
both the topology of the vector B field and in its time evo- 
lution equation, is better appreciated by introducing the 
magnetic vector potential A, defined by B = V x A, so 
that Eq. ([l2]) is automatically satisfied and Eq. ([ll]) takes 
on the following form (in the Coulomb gauge V • A = 0): 



dA 
~dt 



E = 0. 



(13) 



For a given velocity field (the so-called kinematic ap- 
proximation), Eq. (13) may be regarded as a set of 
three-dimensional Hamilton- Jacobi equations, where the 
E components are functions of the spatial first derivatives 
of the A components. Thus, the overall MHD and RMHD 
systems are actually combinations of conservative hyper- 
bolic equations and equations of Hamilton- Jacobi kind. It 
is therefore clear that standard numerical upwind schemes 
developed for Godunov-type hyperbolic sets of equations 
cannot be applied, and proper upwind expressions for the 
magnetic flux functions have to be derived. 

2.2. Characteristic wave speeds in 1-D RMHD 



The one-dimensional case, say d y 



0, is, on the 



V • B = 0, 



(12) 



other hand, almost perfectly equivalent to the hydrody- 
namic case and the overall system can be cast in conserva- 
tive form. In this case Eqs. (]ll]) and ([l2|) yield B x — const 
and Eq. (||) becomes a complete 7x7 system of conserva- 
tion laws, by just adding the [B y , B z ] variables to (||) and 
the fluxes [—E z ,E y ] to ©. However, the 1-D RMHD sys- 
tem, like its MHD counterpart, is not strictly hyperbolic, 
in the sense that two or more eigenvalues may coincide in 
some degenerate cases, depending on the angle between 
the direction of propagation and the local magnetic field. 

The characteristic structure of this system was first 
studied by Anile and Pennisi (1987; see also Anile, 1989), 
who derived the eigenvalues and eigenvectors of the associ- 
ated Jacobian df /du by using the covariant notation. All 
the particular degeneracies were also taken into account. 

However, since in our numerical scheme the detailed 
characteristic eigenstructure is not required, while only 
the two speeds at the local Riemann fan boundaries need 
to computed, here we just report the expressions for the 
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eigenvalues, in the form shown in the appendix of KO. 
These are one entropy wave 



Ao = v x , 

two Alfven waves 



±b x 



i°±b° 



(14) 



(15) 



and four magneto-sonic waves (two fast and two slow 
waves) , that unfortunately do not have a simple analytical 
expression and must be derived from the nonlinear quartic 
equation 

{l-e 2 )(u°X-u x ) 4 

+ (1 - \ 2 )[{c 2 s (b°\ - b x ) 2 - e 2 (u°X - u x ) 2 ] = 0, (16) 
where c 2 = Tp/w is the sound speed squared, b a 



b a /^Wo 



= b a b a 



\b\ 2 /w to t), and e 2 = c 2 s + 



\b\ 2 — c 2 |6| 2 . Note that the ordering of MHD character- 
istic speeds and the various degeneracies are preserved in 
the relativistic case (Anile 1989), although the symmetry 
between each A* couple of waves is lost, due to relativistic 
aberration effects. 

Several numerical algorithms may be employed to solve 
Eq. (jll]): Newton's root finding technique, applied in the 
proper interval for each characteristic speed, Laguerre's 
method for polynomials, involving complex arithmetics, 
or eigenvalues finding routines, based on the associated 
upper Hessenberg matrix. We have tested all these numer- 
ical methods by using the Numerical Recipes (Press et al. 
1986) appropriate routines, which are rtsafe, zroots, and 
hqr, respectively, under a wide range of conditions, includ- 
ing various degenerate cases and ultra-relativistic speeds, 
temperatures or magnetic fields. However, in the code we 
have decided to adopt the analytical approach, described 
for example in Abramowitz & Stegun (1965), which re- 
quires in turn the analytical solution of a cubic and of two 
quadratic algebraic equations. We have found that this 
algorithm gives results comparable to Laguerre's or the 
matrix methods, the most robust and precise ones, and 
it is much less computationally expensive. On the other 
hand, Newton's iterative method, which is the fastest in 
normal conditions, was found to fail in some nearly degen- 
erate cases. 

2.3. Primitive variables 

In order to compute fluxes, the vector of primitive fluid 
variables v = [p, v\ p] T have to be recovered from the 
conservative ones u — [D 1 Qi,E] T , defined in Eq. (||), at 
the beginning of every numerical time step (please notice 
that in the present sub-section E will indicate the total 
energy, which has nothing to share with the electric field 
E). Like for RHD codes, this procedure must be carried 
out by some iterative root-finding routine. In the mag- 
netic case this process is even more difficult, in spite of 



the fact that B can be considered as given (its compo- 
nents are both primitive and conserved variables). In BA 
the full 5x5 system is solved by inverting the u(v) = 
set of nonlinear equations as they stand, providing also all 
the partial derivatives needed. However, we have verified 
that this process is neither efficient nor stable when rela- 
tivistic effects are strong. In Koide et al. (1996) and KO 
the system to solve was reduced down to a couple non- 
linear equations, while here we manage to derive just one 
nonlinear equation to be solved iteratively, with obvious 
improvements both in terms of speed and precision. 

The first step is to use the definitions of u a and b a to 
write the known conservative variables u = [D, Q J , E] T in 
terms of the primitive variables. The vector Q becomes 



Q= (W + B 2 )v — (v ■ B)B, 



(17) 



and by taking the projection along B wc find the impor- 
tant relation S = (Q-B) = W(v-B), where like in Paper I 
we have used W — u>7 2 , w = p + Tip, and Ti = F/(T— 1). 

A 2 x 2 system of nonlinear equations is then derived 
by taking the square of Eq. ( |l7| ) and by using the equation 
for the total energy E: 



(2W + B 2 )B 2 v\ - Q 2 = 0, 



W-p+ is 2 + \b 2 v\ -E = 0, 



(18) 



(19) 



where B 2 v\ = B 2 v 2 - (v ■ B) 2 = B 2 v 2 - S 2 /W 2 . If we 
then use the relations 



P : 



DVl-v^ p=[{l-v 2 )W~p]/T 1 , 



(20) 



it comes out that all quantities appearing in the system 
are written in terms of the two unknowns v 2 (or equiva- 
lently 7) and W. Once these variables are found numer- 
ically, primitive quantities will be easily derived through 
Eqs. and by inverting Eq. (|l7j), that is 



WT&{ Q + W B 



(21) 



In order to bring the system down to just a single 
nonlinear equation (to be solved numerically, for exam- 
ple by Newton's iterative method), we found it useful to 
define B 2 v\ = T 2 /(W + B 2 ) 2 in Eqs. @ and @, where 
T 2 = B 2 Q 2 — S 2 is a new, but given, parameter. Then we 
wr ite Eq. |[|) as a third order algebraic equation for W 
with coefficients that depend on v 2 alone 



1-- 



1-V 



p B T 

w - E+ i + - 



T 2 

(W+B 2 f+ — = 0, (22) 



which can be solved analytically (again, see Abramowitz 
& Stegun 1965). Note that the cubic polynomial on the 
left hand side has a positive local maximum in W = —B 2 . 
Thus, since we know that at least one root must be posi- 
tive, all the three roots of Eq. (^2|) are actually bound to 
be real, and we have verified that the largest one always 
yields the correct result. 



6 



L. Del Zanna et al.: A shock-capturing scheme for relativistic flows - II. 



The function W(£), with £ 



is thus available to- 



gether with its derivative W(£), so the final step is to 
apply Newton's method to find the root of = 0, 

where 



= W 2 £+ {2W + B 2 ) 



and 

^"'(0 = W 2 + 2 WW' 



(W + B 2 ) 2 



T 2 
(W + B 2 ) 



(23) 



(24) 



The numerical routine actually employed in the code is 
rtsafe (Press et al. 1986), which applied with an accuracy 
of 1CT 6 in the range £ min = and £ max = 1 - 1(T 6 (7, nax = 
1000), typically converges in 5-10 iterations for any set of 
conservative variables and magnetic field components that 
actually admits a solution. 

The technique described above appears to be ex- 
tremely efficient and, above all, robust; we therefore rec- 
ommend its use in all RMHD shock-capturing codes, 
whatever the numerical scheme actually employed. 

3. The finite-difference ENO-CT central scheme 

In the present section the third order ENO-CT scheme 
described in LD will be adapted and applied to the 
RMHD equations derived above. The ENO-CT scheme 
employs a finite-difference discretization framework, and 
uses Convex ENO (CENO) reconstruction procedures to 
get high order non-oscillatory point values of primitive 
variables needed to compute numerical flux derivatives. 
The upwind procedures in numerical fluxes are based on 
approximate Riemann solvers of the Lax-Friedrichs type, 
as in other central schemes. In particular, here we adopt 
a flux formula based on two local characteristic speeds, 
rather than just one as in the code discussed in LD. The 
scheme will be here presented in the semi-discrete formal- 
ism, that is time dependency is implicitly assumed for all 
spatially discretized quantities. The evolution equations 
are then integrated in time by applying a third order 
TVD Runge-Kutta algorithm (Shu & Osher 1988), with 
a time-step inversely proportional to the largest (in ab- 
solute value) of the characteristic speeds (the magneto- 
sonic velocities defined in Sect. 2.2) present in the domain 
and subject to the CFL condition. The reader is referred 
to both LD and Paper I for further numerical references 
and comments. In any case, see Shu (1997) for a gen- 
eral overview of ENO schemes, Liu & Osher (1998) for 
the original formulation of the CENO central scheme, and 
Kurganov et al. (2001) for the introduction of two-speed 
averaged Riemann solvers in high order central schemes 
for both hyperbolic and Hamilton- Jacobi equations. 

3.1. Discretization of fluid and magnetic variables 

Given a Cartesian uniform 3-D mesh of N x x N y x N z cells 
(volumes), with sizes Ax, Ay, and Az, let us indicate with 
Pij.k = (xi,yj, Zfc) the cell centers and with P{+i/2,j.k the 



points centered on intercell surfaces (along x in this case). 
Classical Godunov-type schemes are usually formulated 
in a finite-volume (FV) framework, where state variables 
are advanced in time as cell averaged values. By applying 
Gauss' theorem to Eq. (||), the time evolution equation for 
the vector of fluid variables becomes 



df 



Ax 1 

i—l 



at P, 



i,j,k, 



where, for each scalar variable, u 



"1,3 



(25) 



k is the FV cell 



averaged discretization of u(x) at Pi.j t k ■ Flux derivatives 
are given in conservative form as simple two-point differ- 
ences of the numerical fluxes /*. These fluxes are defined 
as intercell surface averages and are stored on surface cen- 
ters, each in the direction corresponding to the i compo- 
nent, with i — 1,2,3. For example, / = f x is located on 
Pi+i/2,j,k points, and the A^ operator, centered in Pi,j,k> 
is defined as 



[A x f]i,j,k — fi+l/2,j,k ~ fi-l/2,j,k- 



(26) 



Similar expressions hold in the other directions. 

A different strategy is needed to reconstruct the mag- 
netic field variables. The induction equation Eq. ([ll]) is in 
curl form, thus the correct procedure for discretization in 
the FV framework is the application of Stokes' theorem. 
The x component gives 



——B r 
dt 



AyE z 

Ay 



A Z Ey 

Az 



at P 



i+l/2,j,k) 



(27) 



where B x is discretized as surface average and located 
at Pi+i/2,j,k intercell points, while for example E z is a 
line average and is located at Pi+x/2,j+i/2,k volume edge 
points (cell corners in 2-D). Similar expressions are de- 
fined for the y and z components. Thanks to the above 
discretization, it is straightforward to prove that the nu- 
merical solenoidal condition will be algebraically satisfied 
at all times (if satisfied at t = 0): 



d_ 

dt 



A X B X 
Ax 



Ay 



A Z B Z 
Az 



= 0, at P, 



i,j,k- 



(28) 



This is the fundamental property of the CT method and 
relies on the definition of the staggered field components 



\B 



x\i+l/2,j,ki [By\i,j + l/2,ki \Pz\i,j,k+l/2 



B- 



(29) 



as primary data. It is important to notice that these data 
contain essential informations not only for the discretiza- 
tion of the corresponding variables at cell centers, but also 
for the definition of longitudinal derivatives at the same 
points, since two values per cell are available for each com- 
ponent. 

A first consequence is that a continuous numerical vec- 
tor potential A can be derived in a unique way by invert- 
ing the discretized form of the V x A = B relation. By 
applying Stokes' theorem once more it is easy to verify 
that these data must be defined as line averages along the 
longitudinal direction and stored at the same locations 
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as the corresponding electric field E components. On the 
other hand, if the A components are used as primary data, 
then the time evolution discretized equations in Eq. ( p7f ) 
must be replaced by 



5* 



-E x at R 



i,j+l/2,fc+l/2) 



(30) 



and similarly for the other components. The field compo- 
nents in Eq. (p9|) are then derived as 



D x — 



A„ 



at R 



i+l/2,j',fe) 



(31) 



and similarly for the y and z components. In the CT 
framework, the two choices are perfectly equivalent, and 
the solenoidal constraint Eq. ( f28| ) is still clearly satisfied 
exactly. 

A second main property is given by the continuity con- 
dition of the face averaged components (i.e. the numerical 
magnetic fluxes) in Eq. (|2^). For example, the numerical 
function 



B x (x) = ^jB x (x)dS, 



(32) 



is a continuous function of x, where S is a cell section nor- 
mal to the x direction. A simple proof is given by integrat- 
ing V • B — on a volume including S and tending to it by 
letting Ax — > (see LD). Therefore, at points of disconti- 
nuity (intercell surfaces) staggered data are well defined as 
point values in the corresponding longitudinal direction, 
having a single-state numerical representation there (this 
fundamental property will be fully appreciated later on for 
upwind calculations) and showing at most discontinuous 
first derivatives along the corresponding coordinate. 

In our ENO-CT scheme, the FV discretization is re- 
placed by a finite-difference (FD) discretization based 
on point values, which is more efficient in the multi- 
dimensional case and allows to use just 1-D interpola- 
tion routines. The primary data actually employed in our 
scheme are thus the point-valued u fluid variables, defined 
at cell centers Pij,k, and the point-valued potential vector 
A components, located at cell edges exactly as in the FV 
approach. The time evolution equations are thus, in the 
FD CT scheme 



i 3 

£« = -E 
dt f-r 


Ax 1 


at Pij^k? 


(33) 


and 








(d/dt)A x = 


-E x , 


a t M,i+i/2,fe+i/2) 




(d/dt)A y = 




at Pi+l/2,j,k+l/2i 


(34) 


(d/d*)A, = 


-E z , 


at -Pj+i/2j+i/2,fe- 





where now the numerical fluid flux functions / are such 
that their volume average approximate the FV fluxes / to 
the given order of accuracy, while electric fields in the FD 
formalism are simply point- valued numerical functions. 



that are integrated by the Runge-Kutta time-stepping al- 
gorithm. The reconstruction procedures and the upwind 

formulae to define numerical fluxes / and E, in the re- 
spective locations, are given in the following sub-sections. 

3.2. Reconstruction procedures 

At this preliminary level of analysis, everything is exact. 
Approximations come into play only in the reconstruction 
procedures, when point values needed for flux computa- 
tions are recovered from primary numerical data to some 
accuracy level. At a second order approximation, the pro- 
cedures are straightforward, since the two discretization 
approaches, FV and FD, coincide. At higher order of ac- 
curacy, on the other hand, specific procedures must be 
defined. 

The first step is to derive magnetic field components 
from vector potential data. This is done in our scheme at 
the beginning of each time-stepping sub-cycle. To third 
order accuracy, we first define 



A x 


= [i 


- 7i^ 2) 


- 7i^ 2) 


Ay 


= [i 


- 7i2?i 2) 


- 7i2?i 2) 


K 


= [i 


- 7i^ 2) 


- 7i^ 2) 



i at -ri+i/2,j+i/2,fej 

where 71 = 1/24 and is the non-oscillatory numerical 
second derivative defined in Paper I. Then the divergence- 
free magnetic field components are derived directly from 
the high order approximation of B = V x A, which gives 



B x = (A y A z )/Ay — (A z Ay)/ Az, at P i+1/ 
B 7 



2,j,k, 



(A Z A X )/Az ~ (A X A Z )/Ax, at P iJ+ i/2 tk , 
(A x A y )/Ax - (A y A x )/Ay, at P iij>k+1 /2- 



(36) 



These new staggered field components clearly satisfy, at 
all times t: 



A X B X Ay By 



A Z B Z 



Aa 



Ay 



0, at Pij t ki 



(37) 



which is the point-value equivalent expression of Eq. (|28|) . 
Thus, in the FD version of the CT scheme, the fundamen- 
tal divergence-free magnetic components are those defined 
in Eq. (|3^), whose divided differences directly give high 
order approximations of the longitudinal derivative. 

Having assured approximated first derivatives satisfy- 
ing exact divergence-free relations, it is now possible to 
reconstruct the corresponding point values in the longi- 
tudinal coordinate, which will be needed for the defini- 
tion of numerical fluxes, in a way to maintain the value of 
the first derivative. Since only one-dimensional operators, 
in turn, are now required, we denote by B the unknown 
point-value numerical data and by B the data derived 
just above, where both sets are located at the same in- 
tercell points. By definition, to third order approximation 
we have 



Eqs. (153) and (34) are thus the time evolution equations [1 — ^{D^\B = B 



(38) 
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For a given set of values B, one has then to invert the 
T>^ operator, which is not based on a fixed stencil of data 
and the resulting matrix appears then to be highly non- 
linear. However, the operator on the left hand side can be 
inverted in an explict way by using the Taylor expansion 

[1 - 7i^ (2) ] _1 =s 1 + [liV (2) ] + [7i 2? (2) ] 2 + ■ • ., so that 
the component B may be approximated by a finite order 
iteration as 



flW = B + YliV^B^ 



BW = B. 



(39) 



It is essential to notice that the approximation order of the 
B point values does not increases with the iteration num- 
ber n, being always of the third order of the base scheme. 
What changes is the residual error in the related solenoidal 
condition: usually n = 5 iterations are enough to assure 
the preservation of the longitudinal derivative value, so 
that Eq. ( |37| ) remains exact within machine accuracy in 
the computation of flux derivatives too, thus avoiding spu- 
rious monopoles terms in the dynamical equations. 

A final interpolation step is needed to define point- 
value magnetic field B components at cell centers Pi,j,k, 
where fluid conservative variables u are stored. Thus, at 
the beginning of each time sub-cycle, point- value primitive 
variables v = [p, v J , p] T can be recovered as described in 
Sect. 2.3. 

In order to use these variables in the definition of nu- 
merical fluxes, a reconstruction step is required along each 
direction, to provide point-value data at intercell points. 
For multi-dimensional calculations and for higher than 
second order schemes, like in our case, the reconstruc- 
tion routines are one-dimensional only in the FD frame- 
work, which is then to be preferred. Moreover, the reason 
for reconstructing the primitive variables is also appar- 
ent: if reconstruction were applied to conservative vari- 
ables, then the time-consuming (in RMHD) algorithm of 
Sect. 2.3 would be needed at intercell points for each di- 
rection. The reconstruction routines employed in our code 
are the Convex ENO procedures described in Paper I, with 
a choice of two slope limiters (MinMod and Monotonized 
Centered) to prevent unwanted oscillations and to pre- 
serve monotonicity. 

Reconstruction procedures at intercell locations give 
two-state, left (L) and right (R), reconstructed variables, 
depending on the stencil used in the definition of the 
(quadratic) interpolation polynomial. However, not all the 
eight variables retain such two-state representation, since 
we know from Eq. (32) that the longitudinal field com- 
ponent along the direction of flux differentiation must 
be continuous at corresponding intercell locations, and 
this property is preserved also for point-value staggered 
B components defined through Eq. (|3^). Therefore, the 
longitudinal field is not reconstructed as the other seven 
Godunov variables, and the single-state value provided by 
Eq. (39) is assumed in flux calculations. Notice that this 
is the crucial point discussed in the introduction: the use 
of the divergence-free field components in numerical fluxes 
permits to avoid the onset of spurious monopoles in the 
computation of the right-hand side of Eq. (|33|). 



Finally, point-value numerical fluxes / defined at in- 
tercell locations (the proper upwind procedures to define 
them will be discussed in the following sub-section) have 
to be further transformed in the corresponding / fluxes 
defined in Eq. ([33]). This step is required to approximate 
flux derivatives to higher than second order. To third order 
accuracy we have, as usual 

7i2? (2) ]/, 



/=[1 



(40) 



whereas this final step is not needed for electric fields, 
which are defined in Eqs. ( |34]) as point-value data. 

3.3. Central-upwind numerical fluxes 

In Lax-Friedrichs central-type schemes, two couples of 
characteristic velocities A± and A^ are first defined at 
intercell locations. These velocities, which are the fast 
magneto-sonic speeds in RMHD, are those at the bound- 
aries of the two local Riemann fans (one fan for each L 
or R state), and are derived by using the procedure de- 
scribed in Sect. 2.2. Then a local averaged Riemann prob- 
lem is solved by using either the two-speed HLL (from 
Harten, Lax, and van Leer) or the single-speed LLF (local 
Lax-Friedrichs) flux formulae: 



/ 



HLL 



pLLF 



where 



„± 



= max{0, ±A ± , ±A R }; a = max{a , a }, 



(41) 



(42) 



(43) 



where basically all the other intermediate states are av- 
eraged out. Notice that when the local Riemann fan is 
symmetric, then a + = aT = a and the two fluxes coin- 
cide, whereas, when both fast magneto-sonic speeds have 
the same sign one of the is zero and the HLL flux be- 
comes a pure upwind flux, either / L or / R . This is why 
the HLL scheme described above was defined in Kurganov 
et al. (2001) as central-upwind. 

The upwind states for the numerical E flux functions 
of the induction equation are defined by the same averaged 
Riemann solver, now to be applied to flux functions having 
a base four-state structure. This different structure arises 
because at cell edges, where the electric fields must be de- 
fined, two surfaces of discontinuity intersect, and modes 
of Riemann fans coming from different directions over- 
lap there. A proper expression that extends Eq. ( [ht| ) and 
( fh2| ) to the induction equation fluxes is derived by tak- 
ing advantage of the analytical and numerical experience 
developed for the Hamilton- Jacobi equations. For ease of 
presentation, only the z component of the electric field, 
E z = —(v x By — v y B x ) (that usually needed in 2-D simula- 
tions), will be treated here, while the x and y components 
are easily obtained by cyclic permutations of the indexes. 

Let us indicate with double upper indexes these four 
states, where the first refers to upwinding along x and the 
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second along y, obtained by reconstructing the required 
primitive variables at the edge point Pi+i/2.j+i/2,k by ap- 
plying a sequence of two independent one-dimensional re- 
construction routines. For each component of the magnetic 
field just one reconstruction in the transverse direction is 
actually required, of course. The proposed HLL and LLF 
upwind formulae for the E z flux function are given respec- 
tively by 

(at + ax)(cty + ay ) 



aj + a x a a y + a 



(B«-B^ x ), (44) 



(45) 



The a x and a^ at -P,:+i/2j+i/2,fc required above should 
be calculated by taking the maximum characteristic speed 
(in absolute value) among the four reconstructed states, 
whereas for sake of efficiency we actually consider the max- 
imum over the two neighboring inter-cell points, where 
these speeds had been already calculated for fluid fluxes. 
As usual, the LLF numerical flux is obtained from the 
HLL one by letting a+ = a~ = a x and a+ = a~ — a y . 

Note that for a pure 2-D case d x A z = —B y and 
d y A z — +B X , thus for a given velocity field the induction 
equation simply becomes the Hamilton- Jacobi equation 
(d/dt)A z = —E z (d x A z ,dyA z ) and our upwind formulae 
correctly match in this case with those given in Kurganov 
et al. (2001), where the HLL central-upwind scheme was 
applied to this kind of equations. Moreover, notice that 
for discontinuity surfaces coincident with one of the inter- 
cell boundaries, that is for 1-D situations, Eqs. ( [5]) and 
( filf ) reduce respectively to Eqs. ( pil"| ) and (|42|), as it should 
be. Thus, for 2-D or 1-D calculations, our 3-D ENO-CT 
scheme automatically treats the magnetic field compo- 
nents which does not require a Hamilton-Jacobi formu- 
lation in the usual Godunov-type approach. 

The fact that the magnetic numerical fluxes for the 
induction equation must use upwind formulae based on 
four-state quantities seems to have been overlooked by 
the other authors, who generally just interpolate the x 
and y 1-D single-state upwind fluxes already calculated 
at intercell points to the cell corner where E z must be 
defined (we specialize here to a 2-D situation). However, 
this procedure is clearly incorrect, because at cell corners 
only B x and B y have a two-state representation, whereas 
v x and v v retain the complete four-state representation. 

4. Numerical results 

The numerical verification of the code is reported here in 
three separate sub-sections: 1-D shock tubes are presented 
in the first, some 2-D test simulations of astrophysical in- 
terest are shown in the second, while the third subsection 



is devoted to more quantitative tests concerning code con- 
vergence on smooth fields in both 1-D and 2-D. 

Since in 1-D RMHD the solenoidal constraint is auto- 
matically satisfied and transverse magnetic field compo- 
nents behave essentially like the other conservative vari- 
ables, the shock tube tests shown here just illustrate the 
ability of the code to handle degenerate cases (where the 
system is no longer strictly hyperbolic due to the coinci- 
dence of two or more eigenvalues) and to separate the var- 
ious Riemann discontinuities or rarefaction waves, which 
are more numerous in the magnetized case (up to seven) 
rather than in the fluid case (just three). On the other 
hand, multidimensional tests truly prove the robustness 
of the code and its accuracy in preserving V • B = in 
time, thus avoiding the onset of spurious forces due to the 
presence of numerical magnetic monopoles. The solenoidal 
constraint is preserved within machine accuracy, like for 
all CT schemes, and therefore the spatial distribution of 
V • B and its evolution in time will not be shown (however, 
see LD for proofs in classical MHD tests). 

For all the simulations that we will show, the scheme 
described above is applied without any additional numer- 
ical viscosity term, which are instead introduced by both 
KO and BA in order to stabilize their Roe- type codes. 
The numerical parameters that may be changed in our 
simulations are just the CFL number (here always c = 
0.5), the reconstruction slope limiter (the Monotonized 
Centered, MC, or the most smearing MinMod, MM, for 
the multidimensional highly relativistic tests shown be- 
low; see Paper I for the precise definition of these lim- 
iters), the order of the reconstruction (third, with CENO 
routines, whenever possible), and the central- type aver- 
aged Riemann solver (we always use the HLL solver). 
Concerning this last point, we have verified that HLL 
and LLF actually give almost identical results in all sim- 
ulations. This is easily explained: in RMHD either the 
sound speed or the Alfven speed are often high, especially 
at shocks where upwinding becomes important, so that 
a + r*> a~ r*> a ~ 1 and the two fluxes tend to coincide. 

4.1. One-dimensional tests 

Shock-tube Riemann problems are not really suited for 
high order shock capturing codes, because oscillations may 
easily appear near discontinuities. This is especially true 
when the reconstruction is not applied to characteristic 
waves, because the various contributions cannot be sin- 
gled out and, for example, it is impossible to steepen 
numerically contact or Alfvenic discontinuities. However, 
we will see here that the base third order CEN03-HLL- 
MC scheme is able to treat this kind of problems reason- 
ably well, usually achieving similar or even better accu- 
racy than characteristics-based second order schemes. In 
Table 1 the parameters for the initial left (L) and right (R) 
states of the proposed Riemann problems are reported (in 
all cases v y = v z = 0, T = 5/3 and t = 0.4). These are 
the same tests as in BA, except the first where T = 2 was 
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0.0 0.2 0.4 O.fl D.B 1.0 0.0 0.2 0.4 OJ O.E 1.0 0.0 0.2 0.4 OJ 0.5 1.0 

Fig. 1. The relativistic analog of Brio & Wu (1988) MHD test problem involving a compound wave. If compared 
to BA, our left moving compound shock and right moving slow shock are better resolved. Here the base scheme 
CEN03-HLL-MC is used, with N = 1600 grid points to compare with BA and Courant number CFL= 0.5. 



Test 


P 


v x 


P 


B x 


By 


B z 


1 L 


1.0 


0.0 


1.0 


0.5 


1.0 


0.0 


1 R 


0.125 


0.0 


0.1 


0.5 


-1.0 


0.0 


2 L 


1.0 


0.0 


30.0 


5.0 


6.0 


6.0 


2 R 


1.0 


0.0 


1.0 


5.0 


0.7 


0.7 


3 L 


1.0 


0.0 


1000.0 


10.0 


7.0 


7.0 


3 R 


1.0 


0.0 


0.1 


10.0 


0.7 


0.7 


4 L 


1.0 


0.999 


0.1 


10.0 


7.0 


7.0 


4 R 


1.0 


-0.999 


0.1 


10.0 


-7.0 


-7.0 



Table 1. Constant left (L) and right (R) states for the 
Riemann problems. 



used. Moreover, for ease of comparison, the same resolu- 
tion used in BA, N = 1600 grid points, is employed. 

The results relative to the first test are shown in Fig. 1. 
This is the relativistic extension of the classic Brio & Wu 
(1988) test, where a compound, or intermediate, shock 
wave is formed. There is still a debate going on about the 
reality of such structures, invariably found by any shock- 
capturing code but not predicted by analytic calculations 
(e.g. Barmin et al. 1996; Myong & Roe 1998). However, it 
is not our intention to contribute to that debate, here we 
just want to show that our third order reconstruction with 
Monotonized Center slope limiter gives better accuracy for 
both the left-going intermediate shock and the right-going 



slow shock, in comparison to the second order scheme of 
BA (which employs a MinMod limiter on magneto-sonic 
shocks and a special steepening algorithm for linearly de- 
generate characteristic variables, i.e. Alfvenic and con- 
tact discontinuities, actually switched off for compound 
waves). As we can see, the total absence of characteristic 
waves decomposition in our code does not prevent at all 
the sharp definition of discontinuities. Moreover, oscilla- 
tions due to high order reconstruction and to the use of 
the most compressive MC limiter, evident in the v x pro- 
file, are kept at a very low level, while, at the same time, 
transitions between constant states and rarefaction waves 
are rather sharp. 

A couple of blast wave examples are shown in Fig. 2, 
again taken from BA, the first with a moderate pressure 
jump and the second with jump as large as 10 4 , produc- 
ing a relativistic flow with a maximum Lorentz factor of 
7 ~ 3.4. Also in these cases our results are basically equiv- 
alent to those in BA, in spite of some small spurious over- 
shoots (more apparent in the p profile in the upper panel 
and in the 7 profile of the second panel), due to the com- 
pressive limiter, and of a rather poorly resolved contact 
discontinuity (in the first test, in the second the density 
peak is far too narrow to recognize it), due to the fact that 
we cannot steepen it artificially because our component- 
wise reconstruction. Again, oscillations are nearly absent 
and rarefaction waves are very well defined. The perfor- 
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iloat wove 2: CENG3-HLL-MC, M = 160G, CFL=0.5 
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Fig. 2. A couple of relativistic, magnetized blast waves. The first (upper panel) has a moderate initial pressure jump 
(p L /p R = 30), whereas the second (bottom panel) has a much stronger jump (p L /p R = 10 4 ), producing a very narrow 
density peak and a Lorentz factor of 7 ~ 3.4. Numerical settings are the same as in Fig. 1. 
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Fig. 3. The relativistic MHD shock reflection problem, with 7 ~ 22.4 flows colliding at x = 0.5. This test is crucial for 
two reasons: the post-shock oscillations, here damped by reconstructing at second order, and the wall heating problem, 
that appears to be quite reduced by the use of HLL. 



mances on this kind of tests mainly depend on the limiter 
choice and on the reconstruction order, so both accuracy 
and numerical problems are similar to those already shown 
in the RHD case. 

Finally, in Fig. 3 we show the fourth test proposed by 
BA, which is the magnetic extension of the shock reflec- 
tion problem of Paper I. To reduce post-shock oscillations, 
more evident in the pressure profile, we have run this test 
at second order, thus as in BA; in spite of this our slow 
shocks are better resolved and the wall heating problem 
produces a lower dip in the density profile at x = 0.5. We 
have also run this test by using highly relativistic flows 
with 7 ~ 224, as in Paper I, and we have not met any 
particular problem. The good performance of our code in 
this last test, in its second order version, is due to the 
use of the HLL solver which is not based on the definition 
of an intermediate state, based on the left and right re- 
constructed quantities, for the definition of characteristic 
speeds, as it is done in usual Roe-type solvers. 

4.2. Multidimensional tests 

For truly multidimensional RMHD tests, analytic solu- 
tions are not available and so the verification of the code 
must be done at a rather qualitative level. Here we will 
present a cylindrical blast explosion, a cylindrical rotat- 
ing disk, both in 2-D Cartesian coordinates with a uni- 



form magnetized medium, and the same astrophysical jet 
of Paper I in cylindrical coordinates, now propagating in a 
magnetized background. The only other 2-D RMHD code 
for wich extensive numerical verification is available in 
the literature is KO, where two blast explosions and a 
Cartesian 2-D jet were tested, in addition to some simula- 
tions of simple 1-D waves and shocks on a 2-D grid which 
will not be repeated here. 

In the first test we use the standard [0, 1] x [0, 1] 
Cartesian grid, here with a resolution of N x = N y — 250 
grid points, and we define an initially static background 
with p = 1.0, p = 0.01 and B x = 4.0. The relativistic flow 
comes out by setting a much higher pressure, p = 10 3 , 
within a circle of radius r = 0.08 placed at the center 
of the domain. Here we use T = 4/3 to reduce plasma 
evacuation at the center. In Fig. 4 we show the situation 
at t — 0.4, when the flow has almost reached the outer 
boundaries. The flow speed reaches its maximum value 
along the x direction, 7, nax ~ 4.35, because the expansion 
of the blast wave is not slowed down by the presence of 
a transverse magnetic field, as it happens along y where 
field lines are squeezed producing the highest magnetic 
pressure. Magnetized cylindrical blast wave are a nice tool 
to investigate the behavior of the plasma, and the robust- 
ness of the code, in a variety of degenerate cases (see KO 
for a detailed description of the various types of shocks 
involved). In this simulation we can see that, despite the 
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Fig. 4. A RMHD 2-D strong explosion with a pressure jump as high as 10 . The resolution is N x — N y = 250 
grid points, and the base multidimensional scheme employing HLL solver and MM limitcr is used. Grayscale levels 
are displayed for density, kinetic pressure, magnetic pressure, and Lorentz factor (together with field lines), where 
5.36 x 10~ 3 < p < 5.79, 0.0 < p < 45.2, 4.32 x 10~ 2 < p m < 72.2, and 1.0 < 7 < 4.35. 



absence of appropriate Riemann solvers handling the de- 
generacies, our code gives smooth and reasonable results 
in all directions. If we compare with the results shown 
by KO, we may see that in spite of the absence of ad- 
ditional artificial resistivity and of the smoothing of the 
initial structure (both included in KO), our results are 
rather smooth. Only low-level noise in the density may 
be seen in the expanding density shell near the diagonals, 
reminiscent of the numerical artifacts found by KO in a 
run without resistivity. These errors are possibly due to 
the use of Cartesian geometry, since numerical errors on 



the independent v x and v y reconstructions are the largest 
precisely along diagonals. 

Another point raised by KO is the possibility of non- 
strict total energy conservation even in CT upwind MHD 
schemes, since magnetic field components are stored and 
evolved at different locations rather than at cell centers 
where fluid variables are defined. However, if the total 
energy, which obviously is a conservative variable, is not 
re-defined in order to prevent unphysical states, it must be 
globally conserved algebraically. We have checked that in 
this 2-D test the total energy is conserved within machine 
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Fig. 5. The relativistic analog of the rotor test, with an initial maximum Lorentz factor of about 10. A high resolution 
(N x — N y = 400 grid points) simulation is shown, with the same numerical settings as in Fig. 4. Grayscale levels 
are displayed for density, kinetic pressure, magnetic pressure, and Lorentz factor (together with field lines), where 
0.35 < p < 8.19, 5.31 x 10~ 3 < p < 3.88, 3.77 x 10~ 4 < p m < 2.43, and 1.0 < 7 < 1.79. 



accuracy, as expected. In our opinion, the results found by 
KO in his set of analogue tests, where the total energy is 
shown to decrease in time (up to a value as large as 3% in 
the intermediate case, see his Fig. 12), are mainly due to 
the presence of a non-consistent treatment of the artificial 
resistivity, which is absent in our code. 

The same numerical parameters, but with a higher res- 
olution (N x = N y = 400), are employed in the second sim- 
ulation, here adapted to the relativistic case from the clas- 
sical MHD one (Balsara & Spicer 1999b; LD; Toth 2000). 
A disk of radius 0.1 with higher density, p = 10, rotating at 



high relativistic speed, u = 9.95 7 max — 10.0, the rotor, 
is embedded in a static background with p = 1.0, p = 1.0 
and B x = 1.0 (r = 5/3). In Fig. 5 the complicated pat- 
tern of shocks and torsional Alfven waves launched by the 
rotor may be seen at the usual output time t — 0.4, when 
the central field lines are rotated of an angle of almost 
90°. This magnetic braking slows down the rotor, whose 
maximum Lorentz at the output time is just 7 = 1.79. 
Note how the initial high density central region has been 
completely swept away: the density has now its minimum 
(p = 0.35) at the center and the material has gone to form 
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Axisymmetric Jet: CEN03-HLL-MM, N,= 400, N,= 160, CFL = 0.5 
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Fig. 6. Magnetized axisymmetric jet simulation in cylin- 
drical coordinates z and r, with N z = 400 and N r = 160 
(20 grid points per inlet radius are employed). The evo- 
lution is shown here for t = 35, with the same settings 
as in Fig. 4 and 5. In the top panel grey-scale shades and 
contours are displayed for —Logio{p), with p m ; n ~ 0.095 
(black) and p max ~ 38.6 (white). In the bottom panel 
Lorentz factor and the magnetic field lines are displayed 
together, where 7 max — 7.14 (black). 

a thin oblong-shaped shell. In spite of the presence of very 
strong shear flows (again, no smoothing is applied to the 
disk boundary in the initial conditions), it appears that 
our central high order HLL solver is good enough both in 
providing high accuracy on smooth waves and in prevent- 
ing numerical oscillations at shocks. The turbulent aspect 
of the high density shell should be due to the nonlinear 
evolution of shear-flow instabilities, since its position co- 
incide with the transition layer where the flow changes its 
direction from tangential to radial. 

Finally, for a truly astrophysical application and as a 
test of the ENO-CT scheme in a non-Cartesian geometry, 
we simulate the propagation of a relativistic axisymmet- 
ric jet in cylindrical coordinates. The initial settings are 
taken the same as in Paper I, for ease of comparison with 
the non-magnetized case. These are a domain [0, 20] along 
z (N z = 400) and [0, 8] along r (N r — 160), corresponding 
to a common resolution of 20 grid points per inlet ra- 
dius, a static background plasma with p — 10.0, p — 0.01, 
B z = 0.1 (r = 5/3), and jet parameters of v z = 0.99 
and p = 0.1, while pressure and magnetic fields are the 
same as in the external medium, corresponding to a den- 
sity ratio n = 1/100 and to a relativistic Mach number 
M = jv/-f CB c s = 18.3 and to a relativistic Alfvenic Mach 
number Ma = "fv/ / y CA CA = 24.3, where c 2 = Tp/wtot and 



c 2 A = B 2 /w to t- Boundary conditions are reflective at the 
axis and extrapolation is assumed across the other bound- 
aries. In the region 0<z<l,0<r<l, the initially 
smoothed jet values are kept constant in time. 

The evolution is shown in Fig. 6 at t = 35, where 
the density logarithm, the Lorentz factor and the mag- 
netic field lines are displayed. Note that the head of the 
jet moves faster than in the non-magnetized case, because 
of the confinement due to the compressed magnetic field 
lines (initial equipartition is assumed, so B 2 = p) that 
also reduces the extension of the cocoon and stabilizes 
Kelvin-Helmoltz instabilities, so nicely defined in Paper I. 
Additional reasons for this latter aspect arc of numeri- 
cal type: the use of MM rather than MC slope limiter and 
the higher numerical viscosity, due to the magnetic contri- 
bution in the fast magneto-sonic speeds, introduces extra 
smoothing of contact discontinuities. 

As we can see from this set of 2-D examples (the 3-D 
case does not present any additional difficulty), our code 
is able to obtain similar results to those shown in KO. 
Like in the fluid case, we have found that the higher order 
of the scheme can compensate the lack of characteristic 
waves decomposition. Even the physical limits that the 
code is capable to cope with look very similar to both 
KO and Koidc ct al. (1996), essentially because of errors 
in reconstruction of multi-dimensional vectors. In fact, 
separate 1-D reconstruction on vector components may 
easily produce unphysical states, for example v 2 may be 
greater than one or the errors on B 2 , needed for fluxes 
and Alfvenic or magneto-sonic speed calculations, may be 
too large leading to states with supcrluminal character- 
istic modes. While the former problem may be cured by 
eliminating the reconstruction in some cases, the latter is 
more difficult to prevent. Another problem, common to 
classical MHD (see Balsara & Spicer 1999a), appears in 
situations of low-beta plasma, where [3 = p/B 2 . When 
the magnetic field is too strong, the pressure is derived 
numerically as a difference of two very large numbers (the 
total energy and the magnetic pressure), so it may even 
become negative. In our code, when the routine described 
in Sect. 2.3 still manages to find a solution for v 2 but then 
negative pressures are found, we reset p to a small value 
(10~ 6 ). The lowest value of the plasma beta that the code 
is able to handle, when relativistic flows are present, ap- 
pears to be around 10~ 3 — 10~ 4 . Typical critical situations 
are strong rarefactions, as in the 2-D blast wave presented 
above. 

4.3. Convergence tests 

The tests shown in the previous subsections were mainly 
devoted to show the robustness of the code on highly rel- 
ativistic flows and shocks, thus also proving the useful 
property of the CEN03 reconstruction algorithm to re- 
duce itself to lower orders near discontinuities, avoiding 
oscillations typical of high-order central- type schemes. In 
the present subsection we check the high-order accuracy 
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CP Alfven wave: CEN03-HLL-MM, CFL = 0.5 
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Fig. 7. Convergence test for the 1-D and 2-D CP Alfven 
wave problem. Relative L\ errors on v z are shown in loga- 
rithmic scale. Note that the code soon achieves third order 
convergence in both cases. The errors in the 2-D case are 
larger because two periods fit on the main diagonal. The 
last point in the 1-D case refers to a slightly larger error 
than expected for perfect third order accuracy, because 
relativistic effects begin to appear on the wave profile and 
we are no longer comparing with the correct solution. 



properties of the CEN03 interpolation routines on smooth 
fields: in cases where discontinuous features are absent, 
these algorithms should actually achieve third order ac- 
curacy. However, the reader might wonder whether the 
overall RMHD scheme, which is based on a rather compli- 
cated sequence of CEN03 reconstruction and derivation 
routines, especially in the multidimensional case, is able 
to preserve global third order accuracy properties in both 
time and space. 

To this porpouse let us study the propagation of rela- 
tivistic circularly polarized (CP) Alfven waves. In the limit 
of small amplitudes the total magnetic field strength is 
preserved in time, the Alfven speed is given by Bq / y/wtot 
and the relation between velocity and magnetic fluctua- 
tions reduce to 5v = ±8By/wtot, similarly to the clas- 
sical MHD case. Define now the various quantities in a 
generic cartesian reference frame (£, 77, Q as p = 1, p = 0.1, 
£?£ = Bo = 1, V£ = 0, and 



Br, 



-Acos(2n£), B ( = -Asin(27rO, 



(46) 



where we have taken A = 0.01. In the 1-D case we sim- 
ply have (£, t], Q — (x,y,z), whereas in the 2-D case we 
consider propagation along the x = y direction, so that 
(£, v, = ((x+y)/V2, {-x+y)/V2, z). In both cases [0, 1] 
intervals and periodic boundary conditions have been as- 
sumed (Bo = V2 in the 2-D case in order to satisfy these 
conditions). 

The convergence can be proved by measuring relative 
errors of a certain quantity, v z in our case, at different 
resolutions, where the error has been evaluated as the L\ 



Li(v z ) 



Taj \vz(xi,yj,t = T) - v z (xi,yj,t = 0)| 
T,n\vz(xi,yj,t = 0)\ 



(47) 



In Fig. 7 the errors are plotted in both 1-D and 2-D 
cases as a function of the number of grid points employed 
N = N x = N y , in logarithmic scale. As expected, third 
order accuracy is achieved already in low resolution runs. 
The base scheme employed is CEN03-HLL-MM, which 
gives the smoothest profiles, more appropriate to wave- 
like features. However, we have also tested the sharper MC 
limiter: third order accuracy is globally preserved, but be- 
haviour of the relative errors as function of the resolution 
appears to be more oscillatory, probably due to artificial 
compression that tends to sharpen somehow even sinu- 
soidal waves (see Fig. 5 of Paper I). 



5. Conclusions 

The shock-capturing 3-D MHD scheme of Londrillo & Del 
Zanna (2000) is applied to the special relativistic case, 
thus extending the code for relativistic gasdynamics de- 
scribed in Del Zanna & Bucciantini (2002), Paper I, to 
the magnetic case. This is the first higher than second or- 
der (third) upwind scheme developed for RMHD, to which 
high resolution Godunov-typc methods have started to 
be applied only very recently. Instead of defining com- 
plicated linearized Riemann solvers, usually based on re- 
constructed characteristic fields, our scheme just uses the 
local fastest characteristic velocities to define a two-speed 
HLL-type Riemann solver. Moreover, reconstruction is ap- 
plied component-wise, thus time-consuming spectral de- 
composition is avoided completely, in the spirit of the 
so-called central schemes. This is of particular impor- 
tance in both MHD and RMHD, since we do not need 
to worry about ubiquitous degenerate cases, usually han- 
dled in Roe-type schemes by adding artificial numerical 
viscosity. 

A main feature of our code is the correct treatment 
of the solenoidal constraint, which is enforced to round- 
off machine errors by extending the constraint transport 
(CT) method, originally developed for the induction equa- 
tion alone, to the overall RMHD system: the flux functions 
are correctly defined by using the staggered magnetic field 
components, thus avoiding the onset of monopoles even 
at discontinuities. It is important to notice that, in or- 
der to obtain such results, V • B must be kept equal to 
zero at cell centers and it must be calculated by using the 
same staggered components which are evolved in time and 
the same discretizations applied to flux derivatives (see 
Toth, 2000, for examples where these properties do not 
apply). Moreover, numerical fluxes based on four-state re- 
constructed quantities are defined for the induction equa- 
tion and here applied to a two-speed central-upwind solver 
for the first time. In our opinion, to date our method is the 
only consistent application of CT to an upwind scheme for 
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mixed systems of hyperbolic and Hamilton- Jacobi equa- 
tions, like MHD and RMHD. 

Particular attention has been also devoted to the nu- 
merical method needed to derive primitive variables from 
the set of conservative ones. The 5x5 system of nonlin- 
ear equation is reduced to just a single equation for the 
square of the velocity. This is then solved via a Newton- 
Raphson iterative root-finding algorithm, and analytical 
expressions are provided for the function whose zeroes are 
looked for and for its first derivative. This procedure is 
extremely efficient and robust, and may be used in all 
RMHD codes, regardless of the numerical scheme em- 
ployed. 

The code is verificated against 1-D shock tube tests 
and 2-D problems, even in non-Cartesian geometries, 
showing accurate results and non-oscillatory profiles. The 
code is very robust within the limits imposed by the intrin- 
sic numerical precision, which for multidimensional rela- 
tivistic flows appear to be 7 ~ 10 — 20 (7 > 200 is reached 
in 1-D calculations) and {3 = p/B 2 ~ 10~ 4 - 10~ 3 . These 
limits seem to be common to all other existing RMHD 
codes, and are essentially due to the fact that physical 
states become undistinguishable in the ultra-relativistic 
regime (e.g. all characteristic speeds collapse onto the 
speed of light), where even very small errors on the re- 
construction produce fluxes that lead to unphysical states. 
Typical situations where the code may fail are strong rar- 
efactions in a strongly magnetized medium. 

Finally, generalized orthogonal curvilinear coordinates 
are defined in the code, and presented in the appendix, 
so our scheme may be easily extended to include General 
Relativity effects with a given metric. 
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Appendix A: Orthogonal curvilinear coordinates 

The equations for special relativistic MHD in a generalized 
orthogonal curvilinear coordinate system (x 1 , x 2 , a; 3 ) are 
obtained by assuming a (covariant) metric tensor of the 
form 



g a p = diag{-l,fr?,/il,/i|}. 



(A.l) 



The first step is to re-define vector and tensor covariant or 
contravariant components as ordinary components, then 
spatial differential operators must be converted in this new 
coordinate system. The set of conservative equations in 
the divergenge form, Eq. (||), becomes 



du 
~dt 



1 



3 



d ( hih 2 h 3 



hih 2 h 3 f—' dx i 



r )+g = o, 



(A.2) 



where u and f l are still those defined in Eq. (||) and (|l0|) , 
respectively. The source term g contains the derivatives of 
the metric elements 
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(A.3) 



where T,j = WtotUiUj — fybj + ptotSij is the stress tensor 



and where we have defined 
1 dhj 



hij = 



hihj dx % 



(A.4) 



Concerning the magnetic evolution equations, since in 
our CT scheme we have assumed A as a primary variable, 
evolved in time by Eq. (|||), the only changes occur in the 
derivation of the magnetic field: 



B, 



\ " 



d 



(h k A k ), 



(A.5) 



which just expresses B = V x A in generalized orthogonal 
coordinates. 

In the code, the various combinations of the metric 
elements are preliminarly calculated and stored on the re- 
quired grids. Thus {hxhzhs) and the six hij terms are 
defined at grid points Pij.k, ^2^3 and its reciprocal are 
defined at the staggered grid Pi+i/2,j,k (similarly for h 3 h\ 
and /11/12), where the corresponding flux and longitudi- 
nal divergence-free magnetic field component need to be 
calculated, and finally the hi elements are stored on the 
same grids where A4 and Ei are defined, that is h\ on 
Pi,j+i/2,k+i/2 a nd so on. Thus, to obtain the derivatives 
in the above expressions, we just need to multiply numer- 
ical fluxes and potential vector components with the cor- 
responding geometrical terms, and then we may proceed 
in the usual way. 
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